Antiferromagnetism of Repulsively Interacting Fermions in a harmonic trap 
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We propose a Real-Space Gutzwiller variational approach and apply it to a system of repulsively 
interacting ultracold fermions with spin ^ trapped in an optical lattice with a harmonic confinement. 
Using the Real-Space Gutzwiller variational approach, we find that in system with balanced spin- 
mixtures on a square lattice, antiferromagnetism either appears in a checkerboard pattern or forms 
a ring and antiferromagnetic order is stable in the regions where the particle density is close to 
one, which is consistent with the recent results obtained by the Real-Space Dynamical Mean-field 
Theory approach. We also investigate the imbalanced case and find that antiferromagnetic order is 
suppressed there. 
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INTRODUCTION 

Ultracold atomic gases have attracted much attention 
[H since the first realization of Bose-Einstein condensa- 
tion 0| . In recent years, ultracold atoms in optical lat- 
tices have stimulated a new wave of studying the many- 
body problems. One can obtain optical lattices by con- 
fining the ultracold atoms in periodic trapping potentials 
created with counter-propagating laser beams [3]. Owing 
to the large degree of control over the optical lattice pa- 
rameters such as the geometry and depth of the potential, 
optical lattices provide an ideal playground for studying 
fundamental condensed-matter physics problems. Many 
remarkable phenomena, like the quantum phase transi- 
tion from a superfluid to a Mott-insulator in a Bose- 
Einstein condensate with repulsive interaction [1] and the 
superexchange interactions with ultracold atoms [H] have 
been observed experimentally in optical lattices. In ad- 
dition, loading ultracold fermions as well as mixtures of 
bosonic and fermionic quantum gases in optical lattices 
has also become a topic of intensive study @, 0, S| ■ 

Although optical lattices have been providing an ideal 
stage for both theoretical and experimental studies 
of fundamental problems in condensed matter physics, 
when compared to true solid state system, defects arise. 
For example, in optical lattices, an additional harmonic 
confinement is always present due to the gaussian profile 
of the laser beams [3j. Although this harmonic confine- 
ment is usually weak and varies slowly (typically around 
10-200 Hz oscillation frequencies) compared to the con- 
finement of the atoms on each lattice site (typically 
around 10-40 kHz), it generally leads to an inhomoge- 
neous environment for the trapped atoms. Therefore, in 
order to make problems more relevant to condensed mat- 
ter systems, investigating how the harmonic confinement 
affects the behavior of atoms trapped in optical lattices 
is important. Motivated by this, we take the ultracold 
fermions with spin \ into consideration and concentrate 
on the magnetic behavior of these particles in such a har- 



monic confinement. 

For simplicity, in this paper we consider the square lat- 
tice with a single orbital per site as a model, which can be 
described by the famous Hubbard Hamiltonian. Hubbard 
model has been studied by various methods such as varia- 
tional Monte-Carlo method 0] and dynamical mean-filed 
theory [Io| . Here we apply the Gutzwiller approximation 
[lH, which was introduced by Gutzwiller along with his 
proposal of Gutzwiller wave function (GWF). It turns 
out that Gutzwiller approximation is exact in the limit of 
infinite dimensions. Extensions to multi-band correlated 
systems using Gutzwiller approximation were carried out 
by J. Biinemann et al. [12|. Meanwhile, Gutzwiller ap- 
proximation wasproved to be equivalent to slave-boson 
theories [HJli, [3] 



on a mean-field level for both one- 



band case [id] and multi-band case 17 



18j. Gutzwiller 
approximation is usually used in homogenous environ- 
ment, here we extends it to inhomogeneous environment 
and address the problem in real space. The organization 
of the paper is as follows: first, we introduce the Hub- 
bard Hamiltonian as well as the Gutzwiller variational 
approach (GVA). Then we show how the harmonic con- 
finement potential and the repulsive interaction affect the 
magnetism of the system in the case of balanced spin- 
mixtures and then we present the results obtained in a 
imbalanced case. Finally, we make some discussions and 
conclusions. 



MODEL AND METHOD 

We apply the Hubbard model for repulsively interact- 
ing fermions in an optical lattice. The Hamiltonian is 
described as 



H — Hq + Hi, 
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where ni„ — 0^0^, and ^(cj^.) are fermionic annihi- 
lation (creation) operators for an atom at the ith site 
with spin a. ty describes the hopping amplitude be- 
tween nearest neighbor sites (ij). If i and j are nearest 
neighbors, iy = t, otherwise, iy =0. U > is the 
repulsive interaction, /v is the chemical potential and 
Vi = \^r\ = Vorf is the external trapping potential, 
in which r, is the distance measured from the center of 
the system. As pointed out in reference 0|, O is usually 
much smaller than the characteristic frequency of the op- 
tical lattice, providing a spatially slowly varying chemical 
potential. 

Many methods, such as Hartree-Fock theory [llj and 
Real-Space Dynamical Mean-Field Theory (R-DMFT) 
approach [2(3], have been used to study the ground state 
of Hubbard model with a confinement potential. Among 
these methods, R-DMFT approach is the most accurate 
and reliable one, because it includes all the local quan- 
tum fluctuation. However, the solution of Anderson im- 
purity model in each iteration step makes it very time- 
consuming. Here we apply the Real-Space Gutzwiller 
variational approach (R-GVA) for this model. We will 
show that the results obtained by R-GVA is consistent 
with those obtained by R-DMFT approach. 

The GVA has been proved to be quite efficient and 
accurate 
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l22l . [23j for the ground state studies of 
many important phenomena in strongly correlated sys- 
tem, i.e. the Mott transition, ferromagnetism and 
superconductivity 24 , 25| . It has also been demonstrated 
[26j that GVA is as accurate as DMFT method for 
the ground state properties, but much computationally 
cheaper, which grants this approach much validity. 

We first give a description of GVA for the ground state 
of correlated electron model systems. There are 2 differ- 
ent spins and each of them could be either empty or oc- 
cupied, thus totally we have 2 2 = 4 number of local con- 
figurations |r) on a single site. Those possible configura- 
tions should not be equally weighted, because electrons 
tend to occupy configurations which have relatively lower 
energy. For this purpose, we could construct projectors 
which can reduce the specified high energy configurations 
|r) on site i 



which fufills, 



m lT = \i,T) (i,T\ 



m iT = 1 



(2) 



(3) 



since all the configurations |T) form a locally complete 
set of basis. 

In Eq.JT]), if the interactions are absent, the ground 
state is exactly given by the Hartree uncorrelated wave 



function (HWF) l^o), which is a single determinant of 
single particle wave functions. However, after turning 
on the interaction terms, the HWF is no-longer a good 
approximation, since it contains many energetically un- 
favorable configurations. In order to describe the ground 
state better, the weights of those unfavorable configu- 
rations should be suppressed. This is the main spirit 
of Gutzwiller wave function (GWF). GWF \^g) is con- 
structed by acting a many-particle projection operator 
on the uncorrelated HWF, 



|*G> = P|*0> 

V =nA:=]lEr^r^r 



(4) 



The projection operator V is used to adjust the weight 
of site configurations through parameters Air (0 < \v < 
1). The GWF falls back to uncorrelated HWF if all 
A^ = 1. On the other hand, if Air = 0, the configu- 
ration r of site % will be totally removed. In this way, 
both the itinerant behavior of uncorrelated wave func- 
tions and localized behavior of atomic configurations can 
be described consistently, and the GWF will give a more 
reasonable physical picture of correlated systems than 
HWF does. 

The evaluation of GWF is a difficult task due to its 
multi-configuration nature. There are lots of efforts in 
the literature, and the most famous one is Gutzwiller 
approximation [ill ]. In this approximation, the inter- 
site correlation effect has been neglected and the physics 
meaning was discussed in [13] and The exact evalu- 
ation of the single-band GWF in one dimension [28| and 
in the limit of infinite dimensions were carried out. 
It turns out that Gutzwiller approximation is exact in 
the latter case. 

The expectation value of Hamiltonian Eq.(l) is: 



(H)g = 



(*g\H\V c 
(*g|*g) 



(* \rHr\y c 

(*o|P 2 |*o) 



(5) 



We note that by choosing Air = y^rr, I*g) is 
normalized under GA. (* g |*g) = IK^ol-Ffl^o) = 



nEi 



(*o|m i; r|*o) = li(Er m *r) = 1. Here 



is the weight of configuration Yi, iriir — (^G\ m ir\^G) 
and m? r = (^ol^irl^o)- In the first equality we sepa- 
rate the average of a projection operator string into the 
product of single site averages. 

The expectation value of kinetic energy is 



(*g|#o|*g) 



where 



(6) 



(7) 
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with -Dp/ r 



<r'| c t |r>|, o 



while for the interaction part of the Hamiltonian 



EE ^<*o|m ir |* > 
i r ll 



(8) 



i r 



Putting Eq.© and Eq.® together, we have the fol- 
lowing equation for the limit of infinite dimensions 

(H)g 

= E t iJ Z ^^ Z J^( C L c ^)a + E £lCT "i,o- + E E K m iT 

i,cr i.V 

= (0\H eff \0) (9) 

In an inhomogeneous systems, the spatial dependence of 
is preserved and the variation is in the 4N S parame- 
ter space, where N s is the number of sites. We adopt the 
following algorithm to minimize (H)q. We begin with 
solving H e ff where the Z-factors are fixed, from which 
we compute the expectation value of the Fermionic oper- 
ators on the ground state. Then the minimization of the 
Gutzwiller variational parameter m is done in the alter- 
nating least squares (ALS) scheme, in which we fix the 
77irs on all but the current site and the problem reduces 
to quadratic optimization which is solved as an eigen- 
problem. Using Eq.JT]), one could compute the Z-factor 
on each site, and then they are plugged into the non- 
interacting model H e ff as parameters. The iteration is 
finished when the difference of Z-factors from two step 
is less then the given precision, say 10~ 6 . In general, 
there is no guarantee that the ALS method will converge 
to the global optimum and the convergence of the itera- 
tion. However, in practice, this does not seem to occur 
as long as one varies the parameters in the Hamiltonian 
adiabatically. 

In the following part, we consider this model on a (24 x 
24) square lattice at half filling (one particle per site) and 
set t as the unit of energy. 



RESULTS AND DISCUSSIONS 

Now we present the numerical results obtained with 
the R-GVA. We focus on the spatial dependence of mag- 
netization and particle distribution at different parame- 
ters. We first consider the balanced situation in which 
the number of spin-t particles is the same to that of the 
spin-| ones. We begin with the discussion on the effect of 
the harmonic confinement V. First we fix the repulsive 
interaction U = 5. The spatial distribution of magneti- 
zation at different strengths of V is shown in Fig. [TJ 
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Figure 1: (Color online) Real-Space magnetization profiles 
for U = 5 on square lattice (24 x 24) at half filling, when (a) 
V = 0.01, (b) V = 0.02, (c) V = 0.03, (d) V = 0.04. AF 
region shrinks as the confining potential increases. 



We find that antiferromagnetic(AF) order exists even 
with the presence of the inhomogeneous harmonic con- 
finement. It is seen clearly from Fig. Q] that how the 
pattern of magnetization evolves as the harmonic con- 
finement V increases. As the confinement potential V 
is enhanced, the antiferromagnetism changes from a uni- 
form checkerboard structure to a ring, where the filling 
is close to 1. 

To make the problem more explicit, we also get the 
particle and spin density profiles along x-direction. In 
Fig. [2] (a) and (b), we present the local density (n)i 
and the absolute value of the staggered magnetism (mi) 
as the function of the distance along y = 12, where 
(m) = (n. iT ) + (mi) and (to,) = ^((n^) - (njj,)). We find 
that in the presence of the confinement potential V, the 
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Figure 2: (Color online) Particle density and the absolute 
value of staggered magnetization as the function of the dis- 
tance along y = 12 for U = 5 at different confining potentials 
at half filling: (a) the density profile; (b) the staggered mag- 
netization. 



antiferromagnetic phase is stable when the local density 
is close to 1 , which is consistent with the result obtained 
in reference ji^l- The results obtained here confirm the 
role that the harmonic confinement plays in affecting the 
antiferromagnetic pattern among the fermions in optical 
lattices. As pointed out in [20j], these results are impor- 
tant for the ongoing attempt to realize antiferromagnetic 
state of fermions with repulsive interactions in periodic 
potentials. 

Next, we concentrate on the effect of the repulsive in- 
teraction U. Experimentally, U could be tuned by the 
Feshbach resonance technique. We first set the confining 
potential V as 0.02. The spatial dependence of magneti- 
zation and the local particle distribution for at different 
strengths of U are presented in Fig. [3] and Fig. [H We 
know that the ground state of ultra-cold fermions loading 
in an optical lattice without trap follows the spin density 
wave (SDW) mean filed prediction at weak coupling. Ap- 
proaching the strong coupling limit, the large repulsive 
interaction drives the system to an AF insulator phase. 
From Fig. [3] and Fig. HI we can see that the confining 
potential V plays a dominant role at weak coupling and 
the SDW state is suppressed, while at strong coupling, 
the repulsive interaction U plays a dominant role and the 
AF order is stable. 

We now investigate the case of imbalanced spin- 
mixtures, i.e. when N-\ ^ N±. The spatial dependence of 
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Figure 3: (Color online) Real-Space magnetization profiles for 
V = 0.02 on square lattice (24 x 24) at half filling, when (a) 
U = 2, (b) U = 3, (c) U = 5, (d) U = 9. The Af region 
expands as U increases. 
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Figure 4: (Color online) Particle density profile along y = 12 
for V = 0.02 at different strengths of repulsive interaction U 
at half filling. The particle density is more and more close to 
1 as U increases. 
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Figure 5: (Color online) Real-Space magnetization profiles for 
V = 0.01 and U = 5 on square lattice (24 x 24) for imbalanced 
spin-mixtures at half filling, when (a) ATj = 288, iVf = 288; 
(b) N l = 286, iV T = 290; (c) N l = 275, iV T = 301; (d) iVj. = 
270, Nf = 306. The AF order decreases as the imbalance is 
enhanced. 




Figure 6: (Color online) Particle density profile as the func- 
tion of the distance along y = x for (7 = 5 and V = 0.01 
for different spin-mixtures at half filling. The particle den- 
sity in the center of the system increases as the imbalance is 
enhanced. 



tering [30| , and by spatial microwave transition and spin- 
changing collisions techniques, which measure the inte- 
grated density profiles along chosen directions (3l| . 



CONCLUSION 

In conclusion, we have developed the fast Real-Space 
Gutzwiller variational approach which is suitable for the 
fast determination of the grounds state phase diagram 
of the inhomogeneous strongly correlated systems. With 
this method, we have studied both balanced and imbal- 
anced case of fermions with spin \ trapped in an optical 
lattice with a harmonic confinement potential. We find 
that the trap potential tends to destroy the AF order in 
the center as well as the edge of the trap, leaving a ring 
of AF region with local density close to 1. The AF order 
is suppressed for imbalanced system. These results are 
meaningful for the ongoing attempt to realize AF in the 
optical lattices. We anticipate that this R-GVA scheme 
could also be applied to other systems, such as a strongly 
interacting Bose- Fermi mixture in a harmonic trap. 



magnetism and the particle density of sublattice at dif- 
ferent strengths of imbalance are presented in Fig. [5] and 
Fig. [6j We find that as the imbalance is enhanced, the 
AF order decreases. In balanced system, antiferromag- 
netism competes with the confining potential V. Upon 
imbalanced spin- mixtures, it follows that an equivalent 
magnetic field is added into the system, therefore the AF 
order is destroyed. 

EXPERIMENTAL SIGNATURES 

Spatial distribution of spin density in a harmonic trap 
predicted in this paper could be detected by Bragg scat- 
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